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I. INTRODUCTION 


A. REASON FOR ANALYSIS 

Many of the satellites currently being launched have highly flexible 
components such as trusses and booms. Long, highly flexible booms and 
trusses are used for gravity gradient stabilization, supporting large solar arrays 
and radio antennas, and ensuring that humans and equipment are adequately 
shielded from radioactive power generators in some proposed interplanetary 
spacecraft. Deploying these trusses and booms and maneuvering the satellite 
for attitude control and station keeping result in large displacements and 
rotations of the tips. This leads to the presence of geometric nonlinearities in 
the structure. Material nonlinearities are not usually present since the booms 
and trusses are not deflected past their yield point and thus the linear stress- 
Strain law applies. Because of the geometric nonlinearities, it is appropriate to 
model flexible space structures by nonlinear equations of motions instead of 
linear ones. Duffing's Equation is one equation used to adequately describe the 
nonlinear behavior due to the elastic effects of these large structures. 

The presence of nonlinearities may result in motions which are impossible 
to predict after long periods of time, yet completely deterministic since the 
equations of motions are known. Slight changes in initial conditions produce 
radically different outputs. These types of responses have been characterized 
as "chaotic." Chaotic vibrations of spacecraft structures could result in control 


problems. 


B. SCOPE OF THESIS 

The purpose of this theses is to demonstrate that chaotic vibrations can 
occur in highly flexible space structures using theoretical and experimental 
analyses. For the theoretical analysis, a flexible sheet of aluminum 7 ft. long, 3 
in. wide, and 0.080 in. thick, was modelled by Duffing's Equations. This process 
is described in Chapter III. For the experimental analysis, the same beam was 
instrumented with strain gauges and accelerometers. The beam was then 
shaken by an MB Dynamics Model PMSOOA Vibration Exciter and the time 
histories of the strains and accelerations were recorded. This process is 
discussed in Chapter IV. Data from both analyses were then analyzed using 
the computer program, "Chaos," developed by CDR Martinus Sarigul-Klijn, 
USN, [Ref. 1]. 

In addition, the purpose of this thesis is to investigate the sensitivity of 
sensors in detecting chaotic vibrations. The two most commonly used types of 
sensors are accelerometers and strain gauges. Measuring accelerations may 
be a more accurate method of recording the time histories of chaotic systems 
since this includes inertia, damping, and stiffness terms of the mathematical 
equations used to model the system. Strain gauges permit measurement of only 
the stiffness terms of the mathematical equations used to model the system. 
After a discussion of the Science of Chaos (Chapter II) the results of the 


theoretical and experimental analyses will be discussed. 


II. THE SCIENCE OF CHAOS 


A. BACKGROUND AND PREVIOUS STUDIES 

A deterministic system whose time history has a sensitive dependence on 
initial conditions is one definition of a chaotic system [Ref. 2: p. 4]. Even though 
the equations of motions which adequately describe a system are well known, 
small changes in the initial conditions produce markedly different outputs. For 
example, although there are many computer programs which can accurately 
predict the weather tomorrow, there are no programs which can accurately 
predict the weather for a day next year. Recognizing that chaos occurs in the 
nonlinear deterministic models raises the hope of understanding the source of 
this random-like behavior and possibly doing something to change it. 

Although Henri Poincare observed chaotic behavior in the trajectories of 
some celestial bodies at the beginning of the twentieth century, the science of 
chaos is a relatively new field. It is one which essentially became possible after 
the birth of the computer which not only allowed numerical integration of 
systems of mathematical equations, but also shorter and shorter time intervals 
at which experimental data could be obtained. In the 1960's Edward Lorenz of 
MIT used a computer to numerically integrate twelve differential equations 
which described a crude model of the weather. In his simulations he found 
dramatic differences in output for subtle changes in input. [Ref. 3: p. 51] Since 
then, many other researchers have used computers to investigate chaotic 
behavior in virtually every field including mathematics, biology, chemistry, 


ecology, economics, physics, medicine, meteorology, and engineering. James 


Gleick presents an interesting and easy-to-read development of the science of 
chaos in Chaos: Making a New Science [Ref. 4]. 

The chaotic behavior of flexible space structures has been studied by 
Moon and Li. In their analysis of a 3.5 meter long 3-D truss, they found that the 
free-play in the joints of the structure resulted in chaotic motion and modal 
frequencies considerably lower than those determined by linear methods. 
Chaotic dynamics in this structure and other space structures might make it 
difficult to design active controls to damp out transient dynamics. [Ref. 5] In 
addition, chaotic dynamics leads to inaccuracies in linear finite element models 
which are currently used to obtain modal analysis during the design of 
spacecraft. The results of the Low Power Atmospheric Compensation 
Experiment (LACE) may lead to improvements of these finite element models. 
[Ref. 6] 

In another study by Moon and Shaw, chaotic vibration were observed in 
an elastic beam with non-linear boundary conditions. For certain values of 
forcing frequency and amplitude, the periodic motion became unstable and 
nonperiodic. The results were typical of a class of mechanical oscillators with 


play or amplitude constraining stops. [Ref. 7] 


B. GEOMETRIC AND TOPOLOGICAL METHODS OF CHAOS 
The purpose of this section is to describe several methods which have been 
developed for detecting chaotic behavior in dynamical systems. The methods 
which are discussed are time series, phase plane, frequency domain, Poincare 
section, Lyapunov exponent, and fractal correlation dimension. The first four 
methods are more qualitative in nature. and answer the question: Is chaos 


present? The last two methods, Lyapunov exponent and fractal correlation 


dimension, are more quantitative in nature. Generally, when determining if 
chaos 1s present in a system, several of the above methods are used since 
applying one method by itself is not always indicative of chaotic motion. 
1. Time Domain Analysis 
The traditional way of observing signals is to view them in the time 
domain. This is simply a record of what is happening to a parameter of the 
system versus time. Sometimes the first indication of chaotic motion 1s that the 
parameter observed exhibits no pattern or periodicity. This test is not always 
conclusive, since the motion observed could have a long period behavior that is 
not readily observable. In addition, some nonlinear systems may appear non 
periodic, but can be broken down into several periodic signals. [Ref. 2: p. 42] 
Figure 2.1 illustrates both the nonperiodic nature of chaotic signals and 
also the sensitive dependence on initial conditions. Thompson and Stewart 


numerically integrated the Duffing's equation 


% + 0.05% + x° =7.5 cost (21) 


to produce the time history in the figure. The ragged appearance persists for as 
long a time as the integrations are carried out. There is never an exact 
repetition of any part of the signal--it is truly nonperiodic. The figure also 
illustrates that for slight changes in initial conditions, the trajectories quickly 


diverge. [Ref. 8: p. 4] 


2. Frequency Domain Analysis 
Because of the problems encountered when trying to detect chaos from - 
the time domain, the frequency domain, or Fourier Spectrum, is primarily used. 
Observing the signal in the frequency domain will instantly reveal if the signal is 
made up of a sum of periodic signals. Chaotic vibrations will exhibit a broad 
spectrum of frequencies in the output when the input isa single frequency 
harmonic motion. Subharmonics and Superharmonics of this frequency will 
appear not as sharp spikes like the excitation frequency but as rounded hills. 
Once again, this test is not always conclusive since multiharmonic outputs may 
result from hidden degrees of freedom of which the observer is unaware. In 
order to detect chaotic vibrations in large degree-of-freedom systems, it is 
necessary to observe changes in the spectrum while one parameter such as 
driving amplitude or frequency is varied [Ref. 2: p. 47]. 
Figures 2.2 and 2.3 show the Fourier spectrum for a periodic signal 


generated by the equation 
X = 3 cos (4x1) (2.2) 


and the chaotic signal of Equation (2.1) respectively. Notice that Figure 2.2 has 

a sharp spike only at the excitation frequency while Figure 2.3 has multiple, 

broadband spikes at frequencies other than the excitation frequency of 0.16 Hz. 
3. Phase Plane 

A dynamical system consists of two parts: the state (essential 

information about a system) and the dynamic (rule that describes how the state 


evolves with time). The motion of the system can be visualized in a state space 


or phase plane. The phase plane is a construct whose coordinates are 
components of the state. The coordinates, therefore, depend on the type of 
system being studied. A mechanical system will usually use position versus 
velocity. A simple example to describe this is the motion of a pendulum. All 
that is needed to determine its motion are position and velocity. The state is 
thus a point on the phase plane. Newton's Laws provide the dynamic, 
expressed mathematically as a differential equation, that describes how the 
state evolves. Without the presence of friction, the pendulum traces out a circle 
or orbit in the phase plane. When friction is present, the orbit spirals in to a 
single point or attractor. [Ref. 3: p. 49] This describes the motion of a simple 
non chaotic system. The orbits of chaotic systems never close or retrace 
themselves. Thus, the trajectory of the orbits will tend to fill up a section of the 
phase plane. Once again, this is not always a conclusive test for chaos. 

The pseudo-phase plane is another technique which has been used 
when only one variable is able to be measured. This technique was discovered 
by scientists at the University of Santa Cruz while analyzing a dripping faucet. 
They constructed a pseudo-phase plane using a measured variable and a time- 
shifted version of the same measured variable. Takens later proved that a 
pseudo-phase plane exhibited the same characteristics as a phase plane. [Ref. 
4: p. 266] The choice of the time shift or embedding time will shrink or 
expand the orbit about the line y = x. Dvorak and Klaschka have proven that 
an embedding time between five and eleven time intervals produces a pseudo- 
phase plane with the least amount of bias [Ref. 9]. Chaotic motions exhibits the 


Same characteristics in the pseudo-phase plane as they do in the phase plane. 


The pseudo-phase plane for the periodic signal appears as a circle as 
expected (Figure 2.4). This signal constantly retraces itself. As seen in Figure 
2.5, the pseudo-phase plane for the chaotic signal of Equation (2.1) never 
retraces itself and fills up a section of the plane. 

4. Poincare Section 

The Poincare section is constructed by passing a two-dimensional 
surface through a three-dimensional phase space and observing where the 
points of the trajectory intersect the plane. This slice can be made through an 
infinite number of locations to reveal the internal structure of the trajectory. 
Periodic signals, such as an undamped pendulum, have one point on the 
Poincare section while signals with multiple harmonics produce a finite number 
of points on the Poincare section. If the Poincare section does not consist of a 
finite set of points or a closed orbit, the motion may be chaotic. In undamped or 
Slightly damped chaotic systems, the Poincare section appears as a cloud of 
unorganized points. In damped chaotic systems, the Poincare section appears 
as an infinite set of highly organized points which look like parallel lines. If after 
several enlargements the structure continues to exist, the system is described 
as having a Strange attractor. This embedding of structure within structure is 
referred to as a Cantor set after the mathematician George Cantor (1845 - 
OS eal. 2 os) 

Figure 2.6 is the Poincare section of the periodic signal. Notice that 
there is only one point in the top right hand corner. The damped, chaotic signal, 
however, appears as a set of points forming parrallel lines as shown in Figure 
2.7. If this figure were enlarged, the same pattern would emerge after each 


enlargement. This is referred to as the self-similar property. 


5. Lyapunov Exponent 

The Lyapunov exponent is a measure of the average rate of 
divergence of trajectories initially separated by an infinitesimal amount. A 
system is chaotic if it has one or more positive Lyapunov exponents, which 
means nearby trajectories are diverging [Ref. 10: p. 323]. Lyapunov exponents 
have proven to be a useful dynamical test for chaotic systems. Since nearby 
trajectories correspond to nearly identical states, exponential divergence 
between the trajectories means that predictability of the system is rapidly lost. 
Systems with subtle, unresolvable initial differences quickly behave differently. 


[Ref. 11: p. 285] The Lyapunov exponent, A, is defined as 


d(t) = d 2" (2.3) 
Or 

ase! d(t) 
where : d is the initial distance between the two trajectories 


d(t) is the distance at a later time, t 

As indicated in Equation (2.4) a positive exponent means that ‘d(t) is greater 
than ‘d’ and the trajectories are diverging. A zero exponent indicates ‘d(t) and 
‘d' are the same and the trajectories coincide. A negative exponent indicates 
‘d(t)' is less than ‘d’ and the trajectories are converging. 

A system with n degrees of freedom has a set of n exponents, called 
the Lyapunov spectrum. This spectrum characterizes the divergence of the 
trajectory inn directions. In order to understand the meaning of the Lyapunov 


spectrum, imagine a sphere of neighboring points in the phase space about a 


central point on the trajectory. As the points evolve with time, this sphere will 
deform. If one axis is diverging exponentially the sphere stretches along that 
direction. If another axis is converging, the sphere shrinks along that direction. 
Algorithms are available to calculate the exponents of the spectrum from the 
measurement of a single observed state by constructing the pseudo-phase 
space. [Ref 12: p. 89] The program "Chaos" uses a method developed by 
Wolf et al (Ref. 11]. 
6. Fractal Correlation Dimension 


The fractal correlation dimension is defined as 
C(r) = r4 (25) 


where: C(r) is the probability that a point on the attractor is within a 
circle, sphere, or hypersphere of radius ‘r' 
ris the radius of the circle, sphere, or hypersphere 
d is the fractal correlation dimension. 
Solving for 'd' in Equation (2.5) yields 


d= lim (TO) (2.6) 


C(r) may be calculated by constructing a sphere at each point on the 


trajectory and counting the number of points in each sphere, such that 


N N 
= res » H (1 - |x; - xj) (2.7) 


1 J 


10 


where: H(s) is the Heaviside function and is equal to one if 's’ is greater 


than zero and equal to zero if ‘s' is less than zero 
| xj - xj |1s the Euclidean distance between points 
N is the total number of points [Ref. 2: p.217] 


The procedure used in the "Chaos" program by Sarigul-Klijn to 


compute the fractal correlation dimension is described below. 


F. 


Ds, 


Start with a point on the attractor and calculate the number of points 
inside a circle of radius, r. 


Calculate the probability C(r) by dividing this number of points by the 
total number of points in the attractor. 


Repeat this for several points along the attractor. 
Compute C(r) for several values of r. 


The slope of log[C(r)] versus log(r) yields d, the fractal correlation 
dimension, for the chosen number of pseudo-phase space dimensions. 
Ref. 1: passi. 


Repeat steps | through 5 for higher pseudo-phase space dimensions to 
obtain an asymptotic estimate of the fractal correlation dimension for the 
dynamical system. 


Nonchaotic systems have integer fractal correlation dimensions such as 


an equilibrium point (zero), limit cycle (one), or toroidal attractor (two). The 


fractal correlation dimension of a chaotic system is non-integer and is 


independent of the dimension of the phase space used for the calculation. [Ref. 


13: p. 101] 
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Figure 2.2 Fourier spectrum of a periodic signal. 
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Figure 2.3 Fourier spectrum of a chaotic signal. 
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Figure 2.4 Pseudo-phase plane of a periodic signal. 
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Figure 2.5 Pseudo-phase plane of a chaotic signal. 
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Figure 2.6 Poincare section of a periodic signal. 
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Figure 2.7 Poincare section of a chaotic signal. 
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Il. THEORETICAL ANALYSIS 


A. DUFFING'S EQUATION 
The mathematical equation which describes the motion of a nonlinear 


elastic beam is represented by the well-known Duffing's equation. 


mV +cv +k,v + k3v? = Acos ot (3.1) 


where: m = mass of the beam 
c = damping coefficient 
kj = linear stiffness coefficient 
k3 = nonlinear stiffness coefficient 
A = amplitude of the excitation 
@® = frequency of the excitation 


v = transverse displacement of the beam 


Once the coefficients of the Duffing’s equation are known it can be integrated 
numerically to solve for the displacement, velocity and acceleration of the tip of 
a cantilever beam. A fourth order Runge-Kutta scheme (see Appendix A) was 
developed to perform the integration. The following methods are used to 
determine the coefficients using the properties of a beam made of 7075-76 


sheet aluminum listed in Table 3.1. 


|e 


B. ASSUMED MODES METHOD 

In order to determine the coefficients to be used in Duffings equation, the 
assumed modes method is used. This method is an extension of the principle 
of virtual displacements and is used to produce a generalized parameter model 
that approximates the flexible behavior of the beam. A virtual displacement is 
an infinitesimal, imaginary change in the displacement of a system as shown in 
Figure 3.1. 

The dotted curve in the figure shows the virtual displacement, dv(x,t), of the 
beam. It must satisfy the same boundary conditions as the actual displacement, 
v(x,t)--the solid curve in Figure 3.1. The virtual displacement is not a function 
of time in the same sense as the actual displacement but is rather a small 
change of configuration relative to the configuration of the beam at the time ‘t™ 
[Ref. 14: p. 33]. 

To create the single degree of freedom model, the deflection of the beam 


may be approximated by 


v(x,t) = wx) Bt) (323 


Similarly, the virtual displacement is given by 
Ov(x,t) = w(x) d0(1) (3.3) 
The shape function, y(x), may be any admissible function which adequately 


describes the behavior of the system and possesses derivatives of an order at 


least equal to those specified by the boundary conditions. The generalized 
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displacement, O(t), is the solution to Duffing’s equation. [Ref. 14: p. 34] The 
shape function which is chosen is the characteristic function for a clamped-free 


beam. 


WG) = ¥2 sin (BS - 7) + exp CF &) + exp (- 5 (1 - §)) (3.4) 


where € is the beam coordinate, nondimensionalized with respect to the length, 
L, and varying between O and 1. [Ref 15: p. 439] 
The boundary conditions for a cantilever beam are 
Bile s(t) = s() 
BC?2: v'(0,t) = 0 


BC3: v (L,t) = 0 
BC4: — v'"(L,t) = 0 


(3.5) 
where (‘) indicates the derivative with respect to the spatial coordinate x. 

The first boundary condition states that there is no displacement at the fixed 
end while the second states that its slope must also vanish here. The third and 
fourth boundary conditions state that there must be no moment and no shear 
forces at the free end. Applying the boundary conditions of Equation (3.5) to 
the shape function of Equation (3.4) gives only a 0.7% error in amplitude at the 
midpoint. [Ref 16: p. 1014] 

The first and second derivatives of the shape function will be necessary in 


determining the coefficients for Duffing's equation. Subsequently, 


"E) =y7E _f).2 _ AL 
y'(§) — 39-7) ag?‘ aml hi 5 (1 - §)) (3.6) 
"E) 7D ren : I sare 
w'"(S) V2" sin Cs ze ? exp ( 5 6) + i exp ( 5 E)) (3.7) 
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The principle of virtual displacements states that the combined virtual work 


of all forces must be zero. More explicitly this is stated as 


dW. - OV TONS 5 + OW. =0 (3.8) 
where OWnc¢ is the virtual work of nonconservative forces, OV is the change in 
Strain energy of the system, OWdamp is the virtual work of damping forces, 
and OWinert is the virtual work of inertia forces [Ref. 14: p. 34]. 
1. Mass Coefficient Determination 
Solving for the virtual work of inertia forces will yield the mass term, m, 


of Duffing's Equation (3.1). Since the inertia force per unit length is distributed 


along the beam, the virtual work of inertia forces is 


L 
OW. ont -| -pA Vv dv dx 
0 


=a | W dé, 8 50 (3.9) 


Using the values for density, area, and length of the beam contained in Table 
3.1 and numerically integrating using MATHCAD (see Appendix B), the work 


done by inertia forces is 


SW inert = 1.0360 x 107! 8 80 (3.10) 


Therefore, the mass coefficient, m, of Duffing's equation is 0.10360 slugs. 
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2. Stiffness Coefficients Determination 
If a beam is modeled as shown in Figure 3.2 with coordinate axes x, y, 
and z and corresponding displacements u, v, and w, it can be shown [Ref. 17: p. 


22} that the strain in the x direction is given by the following equation: 
Exx a HAN ey +o “) ale) 


Assuming that all higher order terms are negligible except for (Ov/dx)2 (small 


Strain, large rotation assumption), the above equation becomes: 


#) 
SoG a. ae (3.12) 


Applying the classical Euler-Bernoulli assumption that cross sections are 
assumed to remain planar and normal to the centroidal axis after deformation 


[Ref. 18: p. 152], it can be seen from Figure 3.3 that 


u = -y tan @, (Gal) 
_ Av 
a ae (3.14) 
Subsequently 
» - Yo (3.15) 
u 
ae ew (3.16) 


Da 


Substituting Equation (3.16) into Equation (3.12) yields 


_ dv , dev” 
Exx = y +o (Gal) 


The expression for the strain energy of a beam undergoing transverse 


deflection is given by [Ref. 19: p. 169} 


L 
vel Oy, Cx Ux On (3.18) 
9) XX “XX 


Because the beam under study is not deflected so the material passes its yield 


point, it behaves as a Hookean material such that 
On = bess (3.19) 


Therefore the strain energy becomes 


vete/ e¢, dx dA (3.20) 
Squaring Equation (3.17) and substituting it into Equation (3.20) yields: 
E 
-E | |yad2vy" _ dvy2d2v 4 Idv ‘as ED) 
v 5 | been’ -y yQy ay + Gy) dx dA 
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If the intersection of the y and z axes lies at the centroid then 


| ydA=L, (6922) 


J vas = (0 (>) 


Therefore, Equation (3.21) becomes: 
_EI] dv" EA | ,dv\4 (3.24) 
V [devia «BAL dota 


Applying variational calculus, Equation (3.24) becomes 


sv =E] | d’vd’6v a, +EA avy’ dby ae (ea) 
dx? dx? : 
which can be written 
oy — ‘| Poy be ae EA | (v')3(8v)' dx (3.26) 


Substituting the derivatives of Equation (3.2) and (3.3) and substituting 
Equations (3.6) and (3.7) into Equation (3.26) yields 


] 


sve Qy")? 0 50 dx + 4 EA | (y')4 o° 50 dx 


0 


(3.27) 


us 


Substituting in the values from Table 3.1 and numerically integrating using 


MATHCAD (See Appendix B), Equation (3.27) becomes 


SV = 1.5433 x 10° 3 80 + 4.5943 x 10! 0° 80 (3.28) 


From Equation (3.28), the linear and nonlinear stiffness coefficients, kjand k3, 
are seen to be 0.015433 Ibs/in and 45.943 Ibs/in3 respectively. 
3. Damping Coefficient Determination 
If the beam is assumed to be modelled by a linear single degree of 


freedom system the damping coefficient is given by 


cq=27¥k)mC (3.29) 
where: C = viscous damping factor 


The values for kj and m were determined previously. In order to 
determine the viscous damping factor, the logarithmic decrement method was 


used [Ref. 14: p. 62]. The logarithmic decrement, 5, is defined by 





8 = In (uP) = BEN (3.30) 
Q eas 
where: Up = amplitude of motion at beginning of cycle 


UQ = amplitude of motion after N cycles 


N = number of cycles 
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The above equation is used to solve for the viscous damping factor by 
displacing the tip of the beam and then measuring the two amplitudes, Up and 
Ug, and counting the number of cycles between the two measurements. These 
measurements were performed on an air table in the Flexible Spacecraft 
Simulator Laboratory of NPS as shown in Figures 3.4 and 3.5. The beam was 
horizontally mounted on the table and an air pad was attached to the tip of the 
beam. The tip was then deflected and the resulting motion achieved essentially 
without friction. The values for Up, Ug, and N appear in Table 3.2. Using 
Equation (3.30), the average viscous damping factor, €, was determined to be 
6.67x10-3 or approximately 0.7%. 


The virtual work of the damping forces is given by 


E 
OW aap = | CqV Ov 
0 


=| yw? dé 3 60 (r3ly) 


Substituting in the value for cq obtained from Equation (3.29) and numerically 


integrating once again yields 


SWaamp = 2.5654 x 10°3 8 59 G2?) 


Therefore, the damping coefficient, c, of Equation (3.1) is 0.0025654 slugs/s. 
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C. THEORETICAL RESULTS 

Appendix A contains the FORTRAN computer code used to numerically 
integrate the Duffing's equation which mathematically describes the motion of 
the flexible beam being studied. Table 3.3 summarizes the results of the 
previous section in which the coefficients of Equation (3.1) were determined. 
Equation (3.1) was then integrated using an amplitude of 3 lbs. (the force 
applied by the shaker) and a frequency of 5.6 Hz. The initial displacement and 
velocity chosen were 1 inch and 0 in/s, respectively. The results were then 
formatted to be read by the computer program "Chaos" developed by M. 
Sarigul Klijn [Ref. 1]. 

Qualitative tests for chaos are satisfied by the time series, Fourier spectrum, 
pseudo-phase plane and Poincare section. Figure 3.6 shows the time history of 
the acceleration of Equation (3.1). It shows no noticeable periodic behavior and 
this is confirmed in the Fourier spectrum of Figure 3.7. The power spectrum is 
made up of many frequencies with broad band characteristics. The pseudo- 
phase plane of Figure 3.8 was obtained using an embedding time of 10. This 
embedding time, consistent with the findings of Dvorak and Klaschka [Ref 9], 
yielded the least amount of bias. Notice that the trajectory frequently crosses 
itself and fills up a section of the phase space. The Poincare section, as shown 
in Figure 3.9, is also characteristic of damped chaotic systems with its cloud of 
unorganized points filling a section of the space. It was obtained by passing a 
plane through the three-dimensional pseudo-phase space at z = 0. 

The quantitative tests for chaos were also performed. Figure 3.10 is a plot 
of the average Lyapunov exponent. This figure is obtained by time averaging 


the many computations of Equation (2.4) and then plotting them versus the 
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number of samples, 4096. The average Lyapunov exponent for this system is 
0.949. This exponent is indicative of chaotic systems. 

Figure 3.11 shows the plot from which the fractal correlation dimension is 
determined. As mentioned previously, the fractal correlation dimension is the 
asymptotic value obtained by plotting the slope of log[C(r)] versus log(r) for 
successively higher pseudo-phase space dimensions. It was determined to be 
approximately 2.4. Moon reports a value of 2.5 for his buckled beam 
experiment. [Ref. 2, p. 232]. The non-integer fractal correlation dimension is 


yet another indication of chaos in dynamic systems. 


Zo 


TABLE 3.1 Beam Properties 
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TABLE 3.2 Logarithmic Decrement Values 
















0.00749 
9375 
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TABLE 3.3 Coefficients of Duffing's Equation 








Figure 3.2 Coordinates of a simple beam 
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Figure 3.3 Kinematics of a simple beam 
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Figure 3.4 View 1 of beam on air table 
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Figure 3.8 View 2 of beam on air table 
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Figure 3.6 Time series of Duffing's 
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IV. EXPERIMENTAL ANALYSIS 


A. EXPERIMENTAL SET-UP 

The aluminum beam with properties listed in Table 3.1 is instrumented with 
Strain gauges and an accelerometer as shown in Figure 4.1. The type of strain 
gauges used are one-quarter inch, 350 ohm, straight geometry gauges from the 
Micro Measurements Group. The accelerometer is a Sensotec model JTF flat 
pack. This accelerometer has a weight of one ounce, a peak range of +5 g's, a 
sensitivity of 10 mV/g, and a usable frequency range of DC to 150 Hz. 

The beam is mounted vertically to an MB Dynamics Model PMS500A 
Vibration Exciter. The outputs from the strain gauges and accelerometer are 
then fed into a bridge balancing box and amplifier. The bridge box is capable of 
handling four inputs, one from the accelerometer and three from the Strain 
gauges. From the bridge box, the signal goes to an analog-to-digital converter. 
From there the acceleration and strain data is stored on a Zenith computer. 
Mechanical stops placed 58 inches up from the base prevent the beam from 


buckling past one and three-fourths inches on either side. (Figure 4.2) 


B. EXPERIMENTAL PROCEDURE 

The accelerometer is placed at location Al in Figure 4.1. The 
accelerometer and the three strain gauges are then connected to the data 
acquisition system. The beam is held in the upright vertical position in order to 
zero Out all the readings from the gauges and the accelerometer using the 
bridge balancing circuit. The beam is then released while simultaneously 


Starting the exciter. The exciter provides an excitation frequency of 5.6 Hz and 


3a 


a displacement amplitude of one-quarter inch peak-to-peak. This procedure is 
then repeated twice, placing the accelerometer at locations A2 and A3. Strain 


gauge measurements at location S1, S2, and S3 are obtained for each run. 


C. EXPERIMENTAL RESULTS 

The data from each of the runs is transferred into a format readable by the 
"Chaos" program. It is then analyzed using the techniques discussed in 
Chapter IJ. These results appear at the end of this chapter. As shown in 


Figures 4.3 through 4.10, the four qualitative tests for chaos seem to be met: 


1. Erratic and unpredictable time histories from both accelerometers and 
Strain gauges. 


2. Multiple broadbanded power spectrum spikes. 
3. Crossing trajectories and filled up pseudo-phase space. 


4. A cloud of unorganized points in the Poincare section characteristic of 
lightly damped systems. 


These results appear similar to the chaotic results obtained by Sarigul-Klijn 
[Ref. 1] in his analysis of helicopter flight test data using higher harmonic 
control. 

The quantitative tests show that the location of the type of sensor chosen 
impacts whether or not chaotic motion is found. Strain gauges only indicate 
chaos if placed near the root, while accelerometers only indicate chaos if 
placed near the tip. Figures 4.11 through 4.14 are plots of the average 
Lyapunov exponents and fractal correlation dimensions. A summary of the 
Lyapunov exponents and fractal correlation dimensions for each sensor appear 
in Tables 4.1 and 4.2. The Lyapunov exponents are all positive for the strain 
gauges. For the accelerometers, the Lyapunov exponent diminishes as the 


accelerometer is placed closer to the root. This indicates that measuring 
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accelerations near the root reduces the chances of obtaining chaotic results. 
This is anticipated since the beam does not experience large deflections in this 
area and can be characterized by linear equations. Likewise, the fractal 
correlation dimensions of the sensors indicate that accelerometers should be 
placed closer to the tip and strain gauges closer to the root. The accelerometer 
nearest the root has an integer fractal correlation dimension, like that of a non- 
chaotic system. The two nearer to the tip have non-integer fractal dimensions, 
characteristic of chaotic systems. Conversely, the two strain gauges near the 
root have non-integer fractal dimensions, while the one closer to the tip has an 


integer fractal dimension. 


a 


TABLE 4.1 Lyapunov exponents 





TABLE 4.2 Fractal Correlation Dimensions 
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Figure 4.12 Average Lyapunov exponents of accelerometers 


(a) Al (b) A2 (c) A3 


=) 


FC.D. 








6 8 10 12 I4 
NUMBER OF DIMENSIONS NUMBER OF DIMENSIONS 
(a) (b) 
0 2 4 6 8 10 12 14 
NUMBER OF DIMENSIONS 
(c) 
Figure 4.13 Fractal correlation dimension of strain gauges 


(ay SP (©) S2 © ss 


52 


FCO. 





0 2 4 6 8 10 12 14 
NUMBER OF DIMENSIONS NUMBER OF DIMENSIONS 
(a) (b) 





0 2 4 6 8 10 il7) 14 
NUMBER OF DIMENSIONS 
(c) 
Figure 4.14 Fractal correlation dimension of accelerometers 


(a) Al (b) A2 (c) A3 


a8 


V. CONCLUSIONS AND SCOPE FOR FUTURE 
RESEARCH 


For the first time geometric and topological methods of chaos were applied 
to data obtained simultaneously by accelerometers and strain gauges. It was 
determined that data from both accelerometers and strain gauges provided 
similar indications to the presence of chaos when analyzed using time series, 
Fourier spectrums, pseudo-phase planes, and Poincare sections. Using these 
techniques, both types of sensors indicated chaos without regard to the location 
of the sensor. However, the Lyapunov exponents and fractal correlation 
dimensions indicated that chaotic motion is only detected by accelerometers if 
they are placed near the tip and that chaotic motion is only detected by strain 
gauges if the are placed near the root. 

The presence of chaos in highly flexible spacecraft structures was also 
investigated. It was determined that the geometric nonlinearities of these 
structures may result in chaotic vibrations and dynamic instabilities. These 
instabilities may then affect control mechanisms of the spacecraft. Further 
analysis could be performed in using the geometric and topological methods of 
chaos to analyze the effectiveness of these mechanisms to adequately control 
the spacecraft and provide improvements into the way the control is performed. 
Chaotic methods may also have an application in updating the linear finite 
elements models of flexible space structures with more accurate nonlinear 


models. 
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TITLE; DUPFFING 
DATE s 07 DEC 1991 
AUTHOR: ROBERT G. VAUGHAN 


SYSTER: VAK 
COMPILER: WASP 


DESCRIPTION: 
This program eolves Duffing’s equation 


mea°2/dt°2(x) + c#d/dt(x) ¢ kiex + k3*x°3 = Acos(w*t) 
by transforming it into two ordinary differential equations 

d/dt(z) = y 

d/at(y) = Acos(wtt) - cty = klex = x3#x73 


and using a fourth order Runge-Kuttsa Scheme to integrate to 
find the displacement, x, and the velocity, y. 


The amplitude and frequency of the forcing function, the ste 
Size for integration, the number of integration steps, and ths 
initial values for displacenent and velocity are entered in 

that order in a data file titled "DUPPFINIT. The data file 
titled "pUFFOUT® econtaing the integratad values for displacement 
and velocity in a format compatible with the program "CHAOS” 

by A. Sarigul-Klijn, 


To alter this program to solve an equation other than Duffing’s 
ubowe, the parameter, NE, must be changed to correspond to the 
nusber of ODE's and the gubroutine "DIFFEQ" must be changed 

to include the proper equations. 


fwo cubroutines are employed. 
DIFFEQ calculates the values for the ODE’s 
RK4 is the Runge-Kutta routine 


As written the program does not use an adaptive etep-size 


The variables are defined as follows: 
NR ~ the number of ODE’s in the systen 
NS ~- the number of steps (itscations) to be performod 
H -— the desired stepeize 
A - amplitude of the forcing funetion 
F —- fraquency of the forcing function 
TIME, T - time 
I, INDEX - integer counters 
X= an arrsy corresponding to the independent variablee, 
Loa. ALLOA, A(Z)=T, etc. 
RDOT - an array corresponding to the ODB’s 


PARRAKARARAKTARETAEARARAEEREAAHHHEEAHHEHKKeEKEHeeHHHAKRHHHAKARHKKANARAANALE SHARES 


PROGRAM DUFFING 


PARAMETER (NE-2) 

IMPLICIT DOUBLE PRECISION(A-H,K,M,O-2) 
CHARACTER*®7 DUMMY 

DIHENSION X(NE) ,XDOT(NE) 


»D 


own 


10 


is 
100 


OPEN(7,FILE=’DUFFINIT’ ,STATUS=’OLD’ ) 
OPEN(8,FILE-’DUFFOUT’ , STATUS’ OLD’ ) 


TIME=0.0D0 
DUMMY=' A1234A’ 


READ(7,*) A,¥F,H,NS 
DO & Iel,NE 
READ(7,*) X(1I) 
CONTINUE 
WRITE(6,6) A,F 
WRITE(8,7) H,NS 
WRITE(8,8) X(1),X(2) 
FORMAT (1X,’FORCING FUNCT AMPLITUDE =~ ’,F4.1,5X,’'FREQ = ‘',F6.3) 


FORMAT (1X,’STEP SIZE = ',F6.4,5X,'NUMBER OF STEPS = ’,15) 
PORMAT (1X,’INIT DISP = ',P4,1,S5X,/INIT VEL = ',F4.1) 


CALL DIFFEQ(X,XDOT,A,F,TIME) 
WRITE(8,100) DUMMY, TIME,X(1),X(2),XDOT(2) 


DO 10 INDEX=1,4 


CALL DIFFEQ(X,XDOT,A,F,TIME) 
CALL RK4(TIME,X,XDOT,H, INDEX) 


CONTINUE 

WRITE(8,100) DUMMY,TIME,X(1),X(2),XDOT(2) 
CONTINUE 
PORMAT(A7, F6.3, E11.4, 2E12.4) 


STOP 
END 


AHKEARARAAAATAKRRAKTAKRHERKAKEEREKEAKAKARERAKAKEKERERAREKEERARERKRKHKAMEKKKKRKAKKE 
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SUBROUTINE RK4(T,X,XDOT,H, INDEX) 


This subroutine performs numerical integration using 4th order 
Runge-Kutta methode. It was developed by I.M. Ross. 


PARAMETER (NE#2) 
IMPLICIT DOUBLE PRECISION(A-H,K,M,O-Z) 
DIMENSION X(NE),XDOT(NE),SAVED(NE), SAVEX(NE) 


GO TO (1,2,3,4), INDEX 


DO 10 I=1,NE 
SAVEX(I)=#X(1I} 
SAVED(I)=XDOT(I) 
X( I} mSAVEX(I)4+0.5D0*H*xXDOT(I) 
T=T+0.5D0*H 
RETURN 


56 


2 DO 20 I#1,NE 
SAVED(1)-SAVED(1I)+2.0D0*xXDOT(I) 


20 X(I)=SAVEX(I)+0.5DO0eH*XDOT(TI) 
RETURN 
3 DO 30 I#1,NE 
SAVED(I)-SAVED(I)+2,0D0*xXDOT(1I) 
20 X(I)—_SAVEX(I)4+H*XDOT(TI) 
T=T+0.5D0*H 
RETURN 
4 DO 40 I=1,NE 
40 X(I)=SAVEX(1I)+5H/6.0D0*(SAVED(I)+xXDOT(I) ) 
RETURN 
END 


cwthkhbkhhkhbhkbkkhhkhkhakenknhhkn kr hk hhh hkhd dade ede hehehehehe kdb hd hdd heheheh hhh hah 


SUBROUTINE DIFFEQ(X,XDOT,A,F,T) 


This subroutine obtains values for the ODE’s for Duffing’s 
. eguation and must be changed for a new system of equations. 


PARAMETER(NE=2 ) 
IMPLICIT DOUBLE PRECISION(A-H,K,M,O-Z) 
DIMENSION X(NE),XDOT(NE) 


PI=-4.0D0*DATAN(1.0D0) 


* 


Be Careful of units! Kl = LBf/in. K3 = LBf/in*3 
* M = slugs C = slugs/s 
M=1.0360D-01 

C=2.5654D-03 

Kij=-1.5433D-02 

K3=4.5943D+01 


* Divide M by 12 in/ft so that XDOT(2) has units in/s~2 
M=M/12.0D0 


XDOT(1)=#X(2) 
XDOT(2)=(A*COS(F*2.0D0*PI*T)-C/12.0D0*X(2)-K1*X(1)-K3*X(1)**3)/M 


RETURN 
END 


3) 


APPENDIX B 
2 
x —T 
1 “Tt: —- (1-x) 
x a 2 2 
m:= J2-sin m~-- - -| +e +e ax 
2 4 
m= 1.388248302957508 
—T —T1 
1 —-x —- (1-x) 
1 x on nm 2 n 2 
k3 := 2°~-COs |T> = — al) =saee + —-@e 
2 é 2 2 


kl = 8. 451750332536703 


58 


10. 


11. 


LIST OF REFERENCES 


Sarigul-Klijn, M., "Application of Chaos Methods to Helicopter Vibration 
Reduction Using Higher Harmonic Control," Ph.D. Dissertation, Naval 
Postgraduate School, Monterey, California, March 1990. 


Moon, F., Chaotic Vibrations, John Wiley and Sons, New York, 1987. 


Crutchfield, J., Farmer, J., Packard, N., and Shaw, R., "Chaos," Scientific 
American, V. 255, December 1986, pp. 46-57. 


Gleick, J., Chaos: Making a New Science, Viking Penguin, New York, 
New York, 1987. 


Moon, F., and Li, G., “Experimental Study of Chaotic Vibrations in a Pin 
Jointed Space Truss Structure," AJAA Journal, V. 28, No. 5, May 1990, 
pp. 915-921. 


Walters, W., "Dynamic Analysis of the Low Power Atmospheric 
Compensation Experiment (LACE) Spacecraft,” Masters Thesis, Naval 
Postgraduate School, Monterey, California, June, 1990. 


Moon, F., and Shaw, S., "Chaotic Vibrations of a Beam with Nonlinear 
Boundary Conditions," /nternational Journal of Nonlinear Mechanics, Vol. 
18, No. 6, 1983, pp. 465-477. 


Thompson, J., and Stewart, H., Nonlinear Dynamics and Chaos, John 
Wiley and Sons, New York, 1986. 


Dvorak, I., and Klaschka, J., “Modification of the Grossberger-Procaccia 
Algorithm for Estimating the Correlation Exponent of Chaotic Systems 
with High Embedding Dimension," Physic Letters A, Vol. 145, No. 5, 16 
Apnil 90, pp. 225-231. 


De Souza-Machado, S., Rollins, R., Jacobs, D., and Hartman, J., "Studying 
Chaotic Systems Using Microcomputer Simulations and Lyapunov 
Exponents," American Journal of Physics, V.58, April, 1990, pp. 321-329. 


Wolf, A., Swift, J., Swinney, H., and Vastano, J., "Determining Lyapunov 
Exponents From A Time Series," Physica D, vol. 16, 1985, pp. 285-317. 


Be 


lish 


16. 


le 


18. 


19. 


Tobochnik, J., and Gould, H., "Quantifying Chaos,” Computers in Physics, 
Nov-Dec, 1989, pp. 86-90. 


Moon, F., and Li, G., "The Fractal Dimension of the Two Well Potential 
Strange Attractor," Physica D, Vol. 17, 1985, pp. 99-108. 


Craig, R., Structural Dynamics, John Wiley and Sons, New York, 1981, pp. 
214-217. 


Dowell, E., "On Asymptotic Approximations to Beam Model Shapes," 
Journal of Applied Mechanics, Vol. 51, June 1984, p. 439. 


Dugundji, J., "Simple Expressions for Higher Vibration Modes of Uniform 
Euler Beams," AJAA Journal, Vol. 26, No. 8, 1988, pp. 1013-1014. 


Oden, J., and Ripperger, E., Mechanics of Elastic Structures, McGraw-Hill, 
New York, 1981. 


Allen, D. and Hailser, W., /ntroduction to Aerospace Structural Analysis, 
John Wiley and Sons, New York, 1985. 


Martin, H., and Carey, G., /ntroduction to Finite Element Analysis, 
McGraw-Hill, New York, 1973. 


60 


INITIAL DISTRUBUTION LIST 


Defense Technical Information Center 
Cameron Station 
Alexandra, VA 22304-6145 


Library, Code 52 
Naval Postgraduate School 
Monterey, CA 93943-5002 


Department Chaiman, Code AA 
Department of Aeronautics and Astronautics 
Naval Postgraduate School 

Monterey, CA 93943-5000 


Department of Aeronautics and Astronautics 
ATTN: Professor Ramesh Kolar, Code AA/Kj 
Naval Postgraduate School 

Monterey, CA 93943-5000 


Department of Aeronautics and Astronautics 
ATTN: Professor I.M. Ross, Code AA/Ro 
Naval Postgraduate School 

Monterey, CA 93943-5000 


Commander M. Sarigul-Klijn, USN 
965 N. Lincoln St. 
Dixon, CA 95620 


Lieutenant Robert G. Vaughan, USN 


11500 Luxmanor Road 
Rockville, MD 20852 


6] 






SHY O19 











Trpasis 


V3525 Vaughan 
al Sensitive u O02 (Seasocs 
for characterizing chaos. 





2 - - 
“ * . =. _- 
- . - . - -- - - 
- 2- - = - ae * -e 
~ -* a om @ * « e 
. . as ” er . 
- = ~ . « « -* . - * , * * . * 
- “ ° PRCR Te oie ee -- o- Peer a te ae oe = . 2 
° . . - . o . “ ’ 
ae oe) asc? - ars woe . ae +82 - ais, Se * -_ . - 
2 « pase ala ~ * ~ . ” - a oe . - - . es @ a . 
PM oan sey Coe iy ood - » wear oF . emee e ' . . - = 
oa rae ‘eg e . ow é . - - . - . aan ne . = e on 
. - « - . oe ~ ~ - o* - . « - 
oe samen o- 6 - ” - eee . ae eo ° . . : ree = 
- -- = « one. 7? » * Paes 7 . Sr eee > =. . . 
aa as - » = ark ee - * s ¢ “= 
- .- . = aoe aces . er ” mot asa apie ‘ 
>» = 6 ._ * 9 . ~- e. - . . - “ . 
--* 7. Paar aac | Se ae eae Se - aa awe * oe , . . 
~ . =< 7s 2 -« . - * ao8 - a * . Py 
._ s- 7? « a zi ee ~~ - el ane ard » . 
~e ~ - . * - a eae oar * . . « + s -« - av 
SC eemater aS e- oe git SAS ee See es - . as ns * - * 
-- « 8 « . - ae * - + . ._ + . ° ° 
=0° <9 @ « a - wert on 8 eo we? . ‘ a © 8 . °- rs . ‘ : 
- - _ * - . * - coer 8 e « =? > > - as 
« nae e . .“e - eae in * - ~ . aoe e « 
eoome © . =. ” = oe ae . . -~* - - *“*e08* - * « bial ° 
“ Pat Week eee ee - . a - . * « = oS ~ 7. - - ace s? . . * 
os Le a aa oes es . . = - Sy te a= * * « 7 fe . a 
. - *° ae -_ = . oe teu. aime: = i. 1a . - °° ~ e .« 1° io a 
. meee PY ee Tet ard « «8 . Pe a ae ee . . @ is 
Sietcts eres Mr te , om - wage a! - ° wrcain o os - « 
«mee -seee ei eerie a owe PCs aye : Se Pere ;, 
> & ‘ .- © - ° jie vee ees «- wre : Ps «< 
e 4 -_e- - ! b ties! J * <? . « ~ ° 
_* * =i ‘ ea AG 
- 1 . 
S - 


ane ave == * 





1 a Bee Beers 





| 
gid 


| 
I 


{| 
{ 
| 
| 


ll 


DUDLEY nin LIBRARY 
3 2768 00018445 1 


ers Bi Cee 
er hee 


et ailed 


6 aye 9 ae aee * 
ead ben tiak reel 
qe getate soeRr OFM # 
mp ee SY 


ee. 


gg lonwteae: 


eOhe ey HS FR 
ae 


eo 
a wena roth 


wot 
eens wewes 





aw OES 
nore 
oteweeour®. 


eect 
uty a 


PEPE yi bak OF ae 
aagaer ere Hike 


i a 
Tis eqrgenner ares 
Cae ne Cee HE WOT ANS 

enequese™ etey 7 
Shake 


Rete 
ohn TKR Ores 
pgesen ate 
ay BERT,” 
ad ALP a oF 


maar ne 
ew ryote 
Fe eee 


cae legere 


eee ee 

= # & 

e « 

- re ee 
. a 





<2 Boe 


= paeete 
+ 0 <8 @.7 
ws 


ga eee 
ne iy J no fae 


vg apy, = OFt Seats bee OS 
< ¢ 


‘ 
Oy vneet oO 


ae ty 


seh Oe 


Gar er Hh Oi Me 


eee § 
wortue, eOe fe 
» yond, Pope MOMs 
ae AHR R Oe © Ge 
Pe bah epee aie 

fame CORTE 


cd 


vou 4 Ss, 
Po 


epeqenente 
hed , 


Pe] 
een ee 
eueenenn 


an 
- 


or) 

~e 

eo te 

wm wr 

ao Bare e F 
“ue n* 





wore hawise 


aan * 


os 
eonet 


oe evar F 
a wet et ees 


rh. ind 
ae 








sar 9 Ah 
wae P We a 

wt Re Dash , .: 
cd eateries 
ee o- - eee Bem, ae ® * 

eee ee eo 

rye Out: 
et ory % 


vans enn * 
ae er 


ee re e oof 
Te gore at ? 
' 9 


iy ad AN PRATT 


ade? od Oe com 


Lag A Re fg ENe 
Segravas © ¥ ? 





16am 8 
teew ota UF a 
ve pcan rere? 

Pratt oat IO 





mr 
woe, Poe 
wee Ont 


met 
waar rere vese 


we Reds Pulls *< 
PS ia 


+e 
wa : 

Pr as ot Soom) 

woreda 8 (Reh? 


PLALY 






are 
argc eter we 


vat ph ery . 
veh Mak ahs 
+e aye 

a oa b 





Wg owed het ¢ C6 A 
Go gh ere ae SS 
wet 


Svat oct iw 
te 
wt 


rene tt 
ee tad 





dase se ener se 
payee te 
hot ie 


wtytes Put 


La 
at beatae tLe 
a tires APs 


eaten 
Per al 
eT ee 
a ba 2 
etal Cat 


agit. 
ae 


+ eeke* # . 
ota heir" = 
yar mers 


Pmezer © 
3°? t,2ar HY 
a ded oe Maer 
oe 2 sen erwin’ * 

w.f Wa" 


at fewe? WP | Poel 
7 +. 


RH Petwsl Fem 


aaa? 





wes 5 
Fempart Eke SWS 
Pruteg  cGedetrte LEN" 
Pq <1 e at ve 


ey Fs 
ate Te 
ape anh S 

. 


met Kee = 


coves by Tee 


“arcs 
he é 


- 
yer at? & 


“« 
ee ee 


ene. a 


i ae, 
rae © 


av 


perv oe Fa 


be = 
. = « 
« e * 
we 
= . « ” 
. +7 
° e - 
oa . 
» » - 
- -_ +> s- . 
- - * 
<4 
s ‘ . road 
- . ot 
= . LY 
vy ~ « 
- . 
« ~*~ 








a = 
,a oe aoe 6 


atete f aot = 
at gts 


+ Ger yest ¢ 
renee : ge Mend od 
mo. te 


tae Ae 


wet. 


hw Pee) 
® 2b ape ree ye? He 

oe wee 
wave 
as 


7-7 
srt coon” 


eT ed 


"4 


eer 
POT Py oe 
pa 


<<» 
oh ater rane ¥ 


pores vepow &* & 
© Gare? oz ett “sy 
‘ hae e+ oe 
ata? Oo 
A epee eee 


wets awh 
feet | ACO = 
- 


eT Sd 


oe 
wat ekegah? 


wet 


a “eter ems 


<onn Oe aa 


yee VES 
aon oe ‘ 


Uke oS 
cerns ¥ 8 


Pit 
> 0, anges dvs 


pity Me 
Ped : 
peer” ade at whe 


me ryte. 
=P oVeen ed ee 
PL tay 


~"¢ 
+ ete te 
ja ee 


priate 
sages see we 


ee aw 


come 
con nee ahs 


Tod Pat Be 

100) “Wy Fee 
Byers ae 8 oe 
otgoad raat 

cate ie ids VERE 
at oMte Reet Ph 

vane 








a 
ae A gh: 
vieewera ates? 
nL ees 

; aet-S 





va, 
se: ehtom VF 
a ods he wed 


ere? 

+ venere 
ete a)" 
sph arge* > 
Lee eee 


ee 
7%? 


ea ar tee 
ore" 


wee aero 


a 





a ® 


qa ot 


nee ee mee 8 





a gs ee er 


| 


- 
- “e - 
ee? . “ oo * 
- = Par 
- Py * 
. 
~ x 
. = - 
« - « 
« 
. “ 
. . 
. . 
- - » = 
ae * * 
* “ =. ” 
ms - o - - ~ 
-- . e 
e . - 
. °? ® our * 
- ° 
. 
a - ~~ . * 
- - 
- « - -—- 
“ . 
2 


ra 
- 


re 


see ta ee 


ee 


wee 
ee eee 


- 


« we -~ 
. . 
- . . 
- - of . * 
7” ow a « - - - 
. = » oa 
w - - . 
. - « 
- - - a 
. * o 2. 
= 
« 
. 
oe 
- 
. <7 





ae 
hy ee ee 


potwsde oe 


a+ 


we 
7 eat eat a wes 
tet 


ware 


wm ou 
sacar oo 


nore 
ape” eet 


matee weer 
on tet eye 2 


gegen ° 
ore e 





eat ee” 


em ees 


yo bs omens 


agiswseeer™ 7" 
wy te ee 


wer 





eatsteee Je 


ewes 


(oa tone, #0 


wees 


Pe dad 


arr 
fase anes 
a hee 


a 
rapt e™ 
+ mee 


ore ce pees 





Saad 
btn v PP Tener eee ~ 8 
xoctpep tere 
eres Sanarvee 
he ay 
pare OL dated 
cath ade one 





agen we HS 


pv care 
Pa 19? 


ogeri ee: 


wa ot aw 
wviae re 
sure? 


are: tatee 
Te bebe 


PP Adah 
= tetere! 
cae ave 


ea ntarnse st. 
fy eeked a Lebel fpeetyms 


=| pee al 
TT Senta c HE gsivecartes 
aeeatedee® 


sat thee. 
rere awest 

meres 

une cots? 






ian eee 


er eo 
rast 


nem 
eos Le 


Ceres 





er ae 


—_ 


acs 


